function [Price, n, V] = euVarTime( r, sigma, S0, K, T, steps, flag, NodeM, D, dnorm )
% euVarTime Summary of this function goes here
% This function uses variable timestepping to calculate PDE of European
% Option prices. Method used is Crank-Nicolson Variable Timestepping

% flag = 0: put (Default), flag = 1: call
% NodeM: Node multiplier, default is 1
% steps: here steps determines the starting dTao = T/steps
% Returned Value: n: time steps

% Argument check
if nargin < 7
    flag = 0; NodeM = 1; D = 1; dnorm = 0.1; end
if nargin < 8
    NodeM = 1; D = 1; dnorm = 0.1; end
if nargin < 9
    D = 1; dnorm = 0.1; end
if nargin < 10
    dnorm = 0.1; end

% Generate S grid
S = [0:0.1*K:0.4*K, 0.45*K:0.05*K:0.8*K, ...
    0.82*K:0.02*K:0.9*K, 0.91*K:0.01*K:1.1*K, ...
    1.12*K:0.02*K:1.2*K, 1.25*K:0.05*K:1.6*K, ...
    1.7*K:0.1*K:2*K, 2.2*K, 2.4*K, 2.8*K, ...
    3.6*K, 5*K, 7.5*K, 10*K];
% Delta Tao
dt = T/steps;
% Number of Asset Prices Grid = Q
Q = length(S);

% Now if we want to increase size of S grid
if NodeM ~= 1
    for i = 2:NodeM
        S_old = S;
        Q = length(S_old);
        S_new = zeros(1,2*Q-1);
        S_new(1:2:end) = S_old;
        S_new(2:2:end-1) = (S_old(1:end-1) + S_old(2:end))/2;
        S = S_new;
        Q = length(S);
    end
end

% Now get the alpha and beta vector
[Alpha, Beta] = alphaBeta(S, r, sigma);

% Now create the sparse matrix M
% This time we use the preallocated indexing sets to create M
% Which is more efficient
i = 2:Q-1;
M = sparse(i,i, (Alpha + Beta + r), Q, Q, 3*Q) + ...
    sparse(i,i-1, -Alpha, Q, Q, 3*Q) + ...
    sparse(i,i+1, -Beta, Q, Q, 3*Q);
M(1,1) = r;

% Option Value at tao = T-T = 0, Default for put option
V = zeros(Q,1);
if flag == 0
    V(Q) = 0;
    V(1:end-1) = max(K - S(1:end-1),0);
else
    V(Q) = S(end);
    V(1:end-1) = max(S(1:end-1) - K, 0);
end
% Now iterate through time steps to solve Vn
% We will use timestep selector to determine delta Tao

% Create sparse identity matrix
I = sparse(1:Q,1:Q,1);

% The C-N Rannacher with variable Timestepping
Tao = dt;   % Current backward time at Tao
dTao = dt;  % initialize dTao
V_old = V;  % initialize Old Option V
n = 1; % Tracker for while loop, total time steps
while Tao < T
    if n <=2  % Fully implicit for first 2 steps
        V_new = (I + M*dTao)\V_old;
    else
        V_new = (I + 0.5*M*dTao)\(I - 0.5*M*dTao)*V_old;
    end
    dTao = varTimeStep(V_new, V_old, dTao, D, dnorm);
    Tao = Tao + dTao;
    V_old = V_new;
    n = n + 1;
end
% Now Tao > T, we do the last step
dTao = T - Tao + dTao;
V = (I + 0.5*M*dTao)\(I - 0.5*M*dTao)*V_old;

% Now interpolate through V to find price based on K
Price = interp1(S,V,K,'nearest');

% Calculate B-S price
[bsCall, bsPut] = blsprice(S0,K,r,T,sigma);
if flag == 0
    BSPrice = bsPut;
else
    BSPrice = bsCall;
end

[delta, gamma] = Greeks(S',V, 2, flag);

end

